home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / RandomArray.py < prev    next >
Text File  |  2006-01-20  |  14KB  |  313 lines
  1. import ranlib
  2. import Numeric
  3. import LinearAlgebra
  4. import sys
  5. import math
  6. from types import *
  7.  
  8. # Extended RandomArray to provide more distributions:
  9. # normal, beta, chi square, F, multivariate normal,
  10. # exponential, binomial, multinomial
  11. # Lee Barford, Dec. 1999.
  12.  
  13. class ArgumentError(Exception):
  14.     pass
  15.  
  16. def seed(x=0,y=0):
  17.     """seed(x, y), set the seed using the integers x, y; 
  18.     Set a random one from clock if  y == 0
  19.     """
  20.     if type (x) != IntType or type (y) != IntType :
  21.         raise ArgumentError, "seed requires integer arguments."
  22.     if y == 0:
  23.         import time
  24.         t = time.time()
  25.         ndigits = int(math.log10(t))
  26.         base = 10**(ndigits/2)
  27.         x = int(t/base)
  28.         y = 1 + int(t%base)
  29.     ranlib.set_seeds(x,y)
  30.  
  31. seed()
  32.  
  33. def get_seed():
  34.     "Return the current seed pair"
  35.     return ranlib.get_seeds()
  36.  
  37. def _build_random_array(fun, args, shape=[]):
  38. # Build an array by applying function fun to
  39. # the arguments in args, creating an array with
  40. # the specified shape.
  41. # Allows an integer shape n as a shorthand for (n,).
  42.     if isinstance(shape, IntType): 
  43.         shape = [shape]
  44.     if len(shape) != 0:
  45.         n = Numeric.multiply.reduce(shape)
  46.         s = apply(fun, args + (n,))
  47.         s.shape = shape
  48.         return s
  49.     else:
  50.         n = 1
  51.         s = apply(fun, args + (n,))
  52.         return s[0]
  53.  
  54. def random(shape=[]):
  55.     "random(n) or random([n, m, ...]) returns array of random numbers"
  56.     return _build_random_array(ranlib.sample, (), shape)
  57.  
  58. def uniform(minimum, maximum, shape=[]):
  59.     """uniform(minimum, maximum, shape=[]) returns array of given shape of random reals 
  60.     in given range"""
  61.     return minimum + (maximum-minimum)*random(shape)
  62.  
  63. def randint(minimum, maximum=None, shape=[]):
  64.     """randint(min, max, shape=[]) = random integers >=min, < max
  65.     If max not given, random integers >= 0, <min"""
  66.     if not isinstance(minimum, IntType):
  67.         raise ArgumentError, "randint requires first argument integer"
  68.     if maximum is None:
  69.         maximum = minimum
  70.         minimum = 0
  71.     if not isinstance(maximum, IntType):
  72.         raise ArgumentError, "randint requires second argument integer"
  73.     a = ((maximum-minimum)* random(shape))
  74.     if isinstance(a, Numeric.ArrayType):
  75.         return minimum + a.astype(Numeric.Int)
  76.     else:
  77.         return minimum + int(a)
  78.      
  79. def random_integers(maximum, minimum=1, shape=[]):
  80.     """random_integers(max, min=1, shape=[]) = random integers in range min-max inclusive"""
  81.     return randint(minimum, maximum+1, shape) 
  82.      
  83. def permutation(n):
  84.     "permutation(n) = a permutation of indices range(n)"
  85.     return Numeric.argsort(random(n))
  86.  
  87. def standard_normal(shape=[]):
  88.     """standard_normal(n) or standard_normal([n, m, ...]) returns array of
  89.            random numbers normally distributed with mean 0 and standard
  90.            deviation 1"""
  91.     return _build_random_array(ranlib.standard_normal, (), shape)
  92.  
  93. def normal(mean, std, shape=[]):
  94.         """normal(mean, std, n) or normal(mean, std, [n, m, ...]) returns
  95.            array of random numbers randomly distributed with specified mean and
  96.            standard deviation"""
  97.         s = standard_normal(shape)
  98.         return s * std + mean
  99.  
  100. def multivariate_normal(mean, cov, shape=[]):
  101.        """multivariate_normal(mean, cov) or multivariate_normal(mean, cov, [m, n, ...])
  102.           returns an array containing multivariate normally distributed random numbers
  103.           with specified mean and covariance.
  104.  
  105.           mean must be a 1 dimensional array. cov must be a square two dimensional
  106.           array with the same number of rows and columns as mean has elements.
  107.  
  108.           The first form returns a single 1-D array containing a multivariate
  109.           normal.
  110.  
  111.           The second form returns an array of shape (m, n, ..., cov.shape[0]).
  112.           In this case, output[i,j,...,:] is a 1-D array containing a multivariate
  113.           normal."""
  114.        # Check preconditions on arguments
  115.        mean = Numeric.array(mean)
  116.        cov = Numeric.array(cov)
  117.        if len(mean.shape) != 1:
  118.               raise ArgumentError, "mean must be 1 dimensional."
  119.        if (len(cov.shape) != 2) or (cov.shape[0] != cov.shape[1]):
  120.               raise ArgumentError, "cov must be 2 dimensional and square."
  121.        if mean.shape[0] != cov.shape[0]:
  122.               raise ArgumentError, "mean and cov must have same length."
  123.        # Compute shape of output
  124.        if isinstance(shape, IntType): shape = [shape]
  125.        final_shape = list(shape[:])
  126.        final_shape.append(mean.shape[0])
  127.        # Create a matrix of independent standard normally distributed random
  128.        # numbers. The matrix has rows with the same length as mean and as
  129.        # many rows are necessary to form a matrix of shape final_shape.
  130.        x = ranlib.standard_normal(Numeric.multiply.reduce(final_shape))
  131.        x.shape = (Numeric.multiply.reduce(final_shape[0:len(final_shape)-1]),
  132.                   mean.shape[0])
  133.        # Transform matrix of standard normals into matrix where each row
  134.        # contains multivariate normals with the desired covariance.
  135.        # Compute A such that matrixmultiply(transpose(A),A) == cov.
  136.        # Then the matrix products of the rows of x and A has the desired
  137.        # covariance. Note that sqrt(s)*v where (u,s,v) is the singular value
  138.        # decomposition of cov is such an A.
  139.        (u,s,v) = LinearAlgebra.singular_value_decomposition(cov)
  140.        x = Numeric.matrixmultiply(x*Numeric.sqrt(s),v)
  141.        # The rows of x now have the correct covariance but mean 0. Add
  142.        # mean to each row. Then each row will have mean mean.
  143.        Numeric.add(mean,x,x)
  144.        x.shape = final_shape
  145.        return x
  146.  
  147. def exponential(mean, shape=[]):
  148.     """exponential(mean, n) or exponential(mean, [n, m, ...]) returns array
  149.       of random numbers exponentially distributed with specified mean"""
  150.    # If U is a random number uniformly distributed on [0,1], then
  151.    #      -ln(U) is exponentially distributed with mean 1, and so
  152.    #      -ln(U)*M is exponentially distributed with mean M.
  153.     x = random(shape)
  154.     Numeric.log(x, x)
  155.     Numeric.subtract(0.0, x, x)
  156.     Numeric.multiply(mean, x, x)
  157.     return x
  158.  
  159. def beta(a, b, shape=[]):
  160.     """beta(a, b) or beta(a, b, [n, m, ...]) returns array of beta distributed random numbers."""
  161.     return _build_random_array(ranlib.beta, (a, b), shape)
  162.  
  163. def gamma(a, r, shape=[]):
  164.     """gamma(a, r) or gamma(a, r, [n, m, ...]) returns array of gamma distributed random numbers."""
  165.     return _build_random_array(ranlib.gamma, (a, r), shape)
  166.  
  167. def F(dfn, dfd, shape=[]):
  168.     """F(dfn, dfd) or F(dfn, dfd, [n, m, ...]) returns array of F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator."""
  169.     return ( chi_square(dfn, shape) / dfn) / ( chi_square(dfd, shape) / dfd)
  170.  
  171. def noncentral_F(dfn, dfd, nconc, shape=[]):
  172.     """noncentral_F(dfn, dfd, nonc) or noncentral_F(dfn, dfd, nonc, [n, m, ...]) returns array of noncentral F distributed random numbers with dfn degrees of freedom in the numerator and dfd degrees of freedom in the denominator, and noncentrality parameter nconc."""
  173.     return ( noncentral_chi_square(dfn, nconc, shape) / dfn ) / ( chi_square(dfd, shape) / dfd )
  174.  
  175. def chi_square(df, shape=[]):
  176.     """chi_square(df) or chi_square(df, [n, m, ...]) returns array of chi squared distributed random numbers with df degrees of freedom."""
  177.     return _build_random_array(ranlib.chisquare, (df,), shape)
  178.  
  179. def noncentral_chi_square(df, nconc, shape=[]):
  180.     """noncentral_chi_square(df, nconc) or chi_square(df, nconc, [n, m, ...]) returns array of noncentral chi squared distributed random numbers with df degrees of freedom and noncentrality parameter."""
  181.     return _build_random_array(ranlib.noncentral_chisquare, (df, nconc), shape)
  182.  
  183. def binomial(trials, p, shape=[]):
  184.     """binomial(trials, p) or binomial(trials, p, [n, m, ...]) returns array of binomially distributed random integers.
  185.  
  186.            trials is the number of trials in the binomial distribution.
  187.            p is the probability of an event in each trial of the binomial distribution."""
  188.     return _build_random_array(ranlib.binomial, (trials, p), shape)
  189.  
  190. def negative_binomial(trials, p, shape=[]):
  191.     """negative_binomial(trials, p) or negative_binomial(trials, p, [n, m, ...]) returns
  192.            array of negative binomially distributed random integers.
  193.  
  194.            trials is the number of trials in the negative binomial distribution.
  195.            p is the probability of an event in each trial of the negative binomial distribution."""
  196.     return _build_random_array(ranlib.negative_binomial, (trials, p), shape)
  197.  
  198. def multinomial(trials, probs, shape=[]):
  199.     """multinomial(trials, probs) or multinomial(trials, probs, [n, m, ...]) returns
  200.            array of multinomial distributed integer vectors.
  201.  
  202.            trials is the number of trials in each multinomial distribution.
  203.            probs is a one dimensional array. There are len(prob)+1 events. 
  204.            prob[i] is the probability of the i-th event, 0<=i<len(prob).
  205.            The probability of event len(prob) is 1.-Numeric.sum(prob).
  206.  
  207.        The first form returns a single 1-D array containing one multinomially
  208.            distributed vector.
  209.  
  210.            The second form returns an array of shape (m, n, ..., len(probs)).
  211.            In this case, output[i,j,...,:] is a 1-D array containing a multinomially
  212.            distributed integer 1-D array."""
  213.         # Check preconditions on arguments
  214.     probs = Numeric.array(probs)
  215.     if len(probs.shape) != 1:
  216.         raise ArgumentError, "probs must be 1 dimensional."
  217.         # Compute shape of output
  218.     if type(shape) == type(0): shape = [shape]
  219.     final_shape = shape[:]
  220.     final_shape.append(probs.shape[0]+1)
  221.     x = ranlib.multinomial(trials, probs.astype(Numeric.Float32), Numeric.multiply.reduce(shape))
  222.         # Change its shape to the desire one
  223.     x.shape = final_shape
  224.     return x
  225.  
  226. def poisson(mean, shape=[]):
  227.     """poisson(mean) or poisson(mean, [n, m, ...]) returns array of poisson
  228.            distributed random integers with specifed mean."""
  229.     return _build_random_array(ranlib.poisson, (mean,), shape)
  230.  
  231.  
  232. def mean_var_test(x, type, mean, var, skew=[]):
  233.     n = len(x) * 1.0
  234.     x_mean = Numeric.sum(x)/n
  235.     x_minus_mean = x - x_mean
  236.     x_var = Numeric.sum(x_minus_mean*x_minus_mean)/(n-1.0)
  237.     print "\nAverage of ", len(x), type
  238.     print "(should be about ", mean, "):", x_mean
  239.     print "Variance of those random numbers (should be about ", var, "):", x_var
  240.     if skew != []:
  241.        x_skew = (Numeric.sum(x_minus_mean*x_minus_mean*x_minus_mean)/9998.)/x_var**(3./2.)
  242.        print "Skewness of those random numbers (should be about ", skew, "):", x_skew
  243.  
  244. def test():
  245.     x, y = get_seed()
  246.     print "Initial seed", x, y
  247.     seed(x, y)
  248.     x1, y1 = get_seed()
  249.     if x1 != x or y1 != y:
  250.         raise SystemExit, "Failed seed test."
  251.     print "First random number is", random()
  252.     print "Average of 10000 random numbers is", Numeric.sum(random(10000))/10000.
  253.     x = random([10,1000])
  254.     if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
  255.         raise SystemExit, "random returned wrong shape"
  256.     x.shape = (10000,)
  257.     print "Average of 100 by 100 random numbers is", Numeric.sum(x)/10000.
  258.     y = uniform(0.5,0.6, (1000,10))
  259.     if len(y.shape) !=2 or y.shape[0] != 1000 or y.shape[1] != 10:
  260.         raise SystemExit, "uniform returned wrong shape"
  261.     y.shape = (10000,)
  262.     if Numeric.minimum.reduce(y) <= 0.5 or Numeric.maximum.reduce(y) >= 0.6:
  263.         raise SystemExit, "uniform returned out of desired range"
  264.     print "randint(1, 10, shape=[50])"
  265.     print randint(1, 10, shape=[50])
  266.     print "permutation(10)", permutation(10)
  267.     print "randint(3,9)", randint(3,9)
  268.     print "random_integers(10, shape=[20])"
  269.     print random_integers(10, shape=[20])
  270.     s = 3.0
  271.     x = normal(2.0, s, [10, 1000])
  272.     if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
  273.         raise SystemExit, "standard_normal returned wrong shape"
  274.     x.shape = (10000,)
  275.     mean_var_test(x, "normally distributed numbers with mean 2 and variance %f"%(s**2,), 2, s**2, 0)
  276.     x = exponential(3, 10000)
  277.     mean_var_test(x, "random numbers exponentially distributed with mean %f"%(s,), s, s**2, 2)
  278.     x = multivariate_normal(Numeric.array([10,20]), Numeric.array(([1,2],[2,4])))
  279.     print "\nA multivariate normal", x
  280.     if x.shape != (2,): raise SystemExit, "multivariate_normal returned wrong shape"
  281.     x = multivariate_normal(Numeric.array([10,20]), Numeric.array([[1,2],[2,4]]), [4,3])
  282.     print "A 4x3x2 array containing multivariate normals"
  283.     print x
  284.     if x.shape != (4,3,2): raise SystemExit, "multivariate_normal returned wrong shape"
  285.     x = multivariate_normal(Numeric.array([-100,0,100]), Numeric.array([[3,2,1],[2,2,1],[1,1,1]]), 10000)
  286.     x_mean = Numeric.sum(x)/10000.
  287.     print "Average of 10000 multivariate normals with mean [-100,0,100]"
  288.     print x_mean
  289.     x_minus_mean = x - x_mean
  290.     print "Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]"
  291.     print Numeric.matrixmultiply(Numeric.transpose(x_minus_mean),x_minus_mean)/9999.
  292.     x = beta(5.0, 10.0, 10000)
  293.     mean_var_test(x, "beta(5.,10.) random numbers", 0.333, 0.014)
  294.     x = gamma(.01, 2., 10000)
  295.     mean_var_test(x, "gamma(.01,2.) random numbers", 2*100, 2*100*100)
  296.     x = chi_square(11., 10000)
  297.     mean_var_test(x, "chi squared random numbers with 11 degrees of freedom", 11, 22, 2*Numeric.sqrt(2./11.))
  298.     x = F(5., 10., 10000)
  299.     mean_var_test(x, "F random numbers with 5 and 10 degrees of freedom", 1.25, 1.35)
  300.     x = poisson(50., 10000)
  301.     mean_var_test(x, "poisson random numbers with mean 50", 50, 50, 0.14)
  302.     print "\nEach element is the result of 16 binomial trials with probability 0.5:"
  303.     print binomial(16, 0.5, 16)
  304.     print "\nEach element is the result of 16 negative binomial trials with probability 0.5:"
  305.     print negative_binomial(16, 0.5, [16,])
  306.     print "\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:"
  307.     x = multinomial(16, [0.1, 0.5, 0.1], 8)
  308.     print x
  309.     print "Mean = ", Numeric.sum(x)/8.
  310.  
  311. if __name__ == '__main__': 
  312.     test()
  313.